Machine Learning–Based Prediction of Disease-Free Survival in Breast Cancer Incorporating Clinical and Treatment Information: Retrospective Cohort Study With External Validation

Introduction

Breast cancer is the most prevalent malignancy worldwide, with the latest data indicating approximately 2.3 million new breast cancer diagnoses annually, making it a leading cause of cancer-related mortality among women [1]. Although the overall survival rate of patients with breast cancer has significantly improved due to the widespread adoption of early screening technologies and advancements in comprehensive treatment plans, including surgery, chemotherapy, radiotherapy, and endocrine therapy, postoperative disease-free survival (DFS) remains a severe challenge for long-term survival [2]. It has been reported that approximately 20% to 30% of patients with breast cancer experience local DFS or distant metastasis within 5 to 10 years after diagnosis, which severely compromises patient survival prognosis and quality of life and imposes a substantial burden on health care systems and society [3]. Therefore, the accurate assessment of DFS risk in patients with breast cancer has important clinical implications. It can provide an intensified therapy protocol for high-risk patients and guide the formulation of individualized adjuvant treatment plans and follow-up strategies.

Currently, the clinical assessment of breast cancer DFS prognosis traditionally relies on the tumor, node, metastasis (TNM) staging system established by the American Joint Committee on Cancer [4]. However, breast cancer exhibits high heterogeneity, implying that patients with similar clinical stages and pathological characteristics can have markedly different clinical outcomes [5]. Although these factors have been incorporated into classic prognostic models such as the Nottingham Prognostic Index, the predictive performance is limited by the heterogeneity of breast cancer and the confounding relationships among multiple prognostic variables [6]. Furthermore, with the development of high-throughput technologies and the availability of multiomics data such as gene expression, DNA methylation, and proteomics, traditional statistical models (eg, Cox proportional hazards regression) demonstrate certain limitations in analyzing complex, nonlinear interactions within clinical data [7].

In recent years, machine learning (ML) algorithms have emerged as powerful tools for processing high-dimensional, nonlinear, and heterogeneous data, offering novel approaches to enhance the accuracy and robustness of DFS prediction models [8]. Unlike traditional statistical methods that require explicit assumptions about the data distribution and variable relationships, ML algorithms can automatically learn complex patterns and latent associations from large-scale datasets through iterative training and optimization [9]. Previous studies have shown that ML models integrating clinicopathological features, multiomics data, or imaging data outperform traditional prognostic models in predicting breast cancer DFS [10]. Additionally, ML models can provide interpretable results through feature importance analysis, aiding in the identification of novel biomarkers and offering new prognostic indicators for clinical practice.

Despite the advantages and general applicability of ML models over traditional prediction methods, their clinical application can still lead to biased results due to inadequate dataset size and training data [11]. Therefore, it is essential to develop robust and interpretable ML models for breast cancer DFS prediction using large-scale sample sets, followed by rigorous external validation and clinical correlation assessment. This study aimed to develop a DFS prediction model that integrates clinicopathological features from The Cancer Genome Atlas (TCGA) and an independent validation cohort to predict breast cancer DFS risk and evaluate its performance relative to traditional prognostic models and existing ML models.

Methods

Data Source and Patients Selection

A total of 1869 female patients diagnosed with breast cancer were retrospectively screened from January 2012 to December 2020 at the Affiliated People’s Hospital of Jiangsu University. Patients were eligible if they met the following criteria: age greater than 18 years, newly diagnosed T1-4 N0-3 M0 breast cancer confirmed by pathological diagnosis, and receipt of surgical treatment. Patients who were lost to follow-up, aged more than 75 years, had distant metastasis at initial diagnosis, or had an uncertain diagnosis were excluded. Finally, 1231 patients were included in the institutional cohort.

The 1231 patients in the institutional cohort were randomly divided into a training cohort (n=861) and an internal validation cohort (n=370) at a ratio of 7:3. The training cohort was used for feature selection, model development, and hyperparameter optimization, whereas the internal validation cohort was used to evaluate model performance in previously unseen patients. The difference in sample size between the training and internal validation cohorts was a direct consequence of the prespecified 7:3 random partitioning strategy.

For external validation, data from 1097 patients were independently obtained from TCGA database, including clinical information and RNA sequencing profiles. Applying the same eligibility criteria, 304 patients were included in the external validation cohort. The TCGA cohort was completely independent of the institutional cohort and was not used for feature selection, model development, or hyperparameter optimization.

Missing data in the TCGA cohort were handled using multiple imputation with the “mice” package in R [12]. For the continuous variable lymph node ratio (LNR; 27/304 missing), predictive mean matching was applied. Binary variables, including menopausal status (28/304 missing), chemotherapy, endocrine therapy, and targeted therapy (32/304 missing), radiotherapy (19/304 missing), and estrogen receptor (ER) status (1/304 missing), were imputed using logistic regression.

In the institutional cohort, Ki-67 expression was assessed by immunohistochemistry and categorized as low (<20%) or high (≥20%), based on the 2013 St. Gallen International Expert Consensus and the 2018 guidelines [13,14]. This 20% cutoff was selected to reflect contemporary clinical practice and consensus-based assessment during the study period. Because standardized Ki-67 and P53 immunohistochemical measurements were not available in the TCGA-BRCA (The Cancer Genome Atlas Breast Invasive Carcinoma) cohort, RNA sequencing data were used as transcriptomic surrogate measures, with MKI67 and TP53 expression corresponding to Ki-67 and P53, respectively. For external validation, the MKI67 and TP53 expression values were dichotomized using thresholds derived to approximate the corresponding distribution in the institutional cohort based on the 25th percentile (lower quartile) of expression. These transcriptomic measures were used as validation surrogates and were not considered directly equivalent to immunohistochemistry-derived Ki-67 and P53 measurements in the institutional cohort. DFS was defined as the time from the date of surgery to the first occurrence of disease recurrence, metastasis, or death from any cause. Patients without any of these events were censored at the date of their last available follow-up.

The candidate predictors included demographic, pathological, perioperative, and postoperative treatment variables. Demographic and pathological variables included age, menopausal status, LNR, tumor laterality, pT stage, pN stage, ER status, progesterone receptor status, and human epidermal growth factor receptor 2 (HER2) status. Surgery type was considered a perioperative clinical variable. Chemotherapy, endocrine therapy, targeted therapy, and radiotherapy represented postoperative treatments actually received by the patients rather than treatment plans available at the time of surgery. These treatment variables were incorporated to develop a postoperative, treatment-informed DFS prediction model. Their predictive importance was interpreted as prognostic information contained in the observed postoperative treatment history and not as evidence of a causal treatment effect.

Model Building and Evaluation

The 1231 patients with breast cancer were randomly divided into a training cohort and an internal validation cohort at a ratio of 7:3. In addition, 304 patients with breast cancer screened from the TCGA database were included as an external validation cohort. Five ML methods, including random survival forest, survival support vector machine, CoxBoost, light gradient boosting machine (LightGBM), and extreme gradient boosting (XGBoost) with an accelerated failure time (AFT) objective, together with Cox regression, were used to develop prediction models. Harrell’s C-index, the area under the curve (AUC) at 3, 5, and 8 years, and the integrated Brier score were used to compare the 6 models. Calibration curves at 3, 5, and 8 years and decision curve analysis were used to evaluate model calibration and clinical utility. In the calibration curves, the number of bins was adaptively adjusted based on the sample size of each cohort to ensure the stability of the Kaplan-Meier estimates within each group [15]. Bootstrap resampling method (2000 iterations) was applied during model development and validation to enhance the reliability of the study results [16].

Model Interpretability and Risk Modeling

The Shapley additive explanations (SHAP) method is an advanced ML tool used to interpret the outputs of predictive models. SHAP analysis was performed on the XGBoost (AFT) model to enhance model transparency and facilitate clinical interpretability. The global importance of each variable was quantified using mean absolute SHAP values, and the variables were subsequently ranked. The SHAP force plot was used to interpret individual model predictions. It decomposes the predicted output into contributions from each feature, starting from the base value. Positive and negative SHAP values indicate features that increase or decrease the predicted outcome, respectively. This plot intuitively shows the magnitude and direction of each feature’s effect on individual prediction results.

Variables with a mean absolute SHAP value contribution of 5% were considered to have a nonnegligible contribution to model prediction and were thus selected as candidate predictors. Subsequently, these candidate factors were used to construct a composite DFS risk score model. The risk score for each patient was calculated as a weighted linear combination of the selected variables and their corresponding regression coefficients. The performance of the risk score model was validated using XGBoost (AFT), and the model was applied for risk stratification and further prognostic assessment. The study was reported in accordance with the guidelines (Checklist 1).

${\displaystyle {\displaystyle {\displaystyle \mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\sum _{i=1}^N({\beta }_i\times X_i),}}}$

where N is the number of prognostic variables, X_i represents the ith variable, and β_i is the regression coefficient derived from the Cox regression model.

Statistical Analysis

The comparison of clinicopathological characteristics between the training and internal validation cohorts was conducted using the Mann-Whitney U test, whereas the comparison of clinicopathological characteristics among the training, internal validation, and external validation cohorts was conducted using the Kruskal-Wallis H test. Given that all variables included in this study were categorical, intervariable correlations were evaluated using Cramér V coefficient to ensure methodological appropriateness and robustness [17].

ML model development was performed using Python (version 3.13.9), incorporating algorithms including XGBoost (AFT), LightGBM, random survival forest, and fast survival support vector machine, implemented via the scikit-survival framework. Detailed hyperparameter settings for each model are provided in Table S1 in Multimedia Appendix 1. All remaining statistical analyses and visualizations were conducted using R (version 4.5.2; R Foundation for Statistical Computing), primarily relying on the survival, tidyverse, ggplot2, tableone, and mclust packages. A 2-tailed P<.05 was considered statistically significant.

Ethical Considerations

This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. All data were anonymized, and the requirement for informed consent was waived due to the retrospective nature of the study. The study was reviewed and approved by the Ethics Committee of the Affiliated People’s Hospital of Jiangsu University (approval ID: [2026] KYW001-27) on January 23, 2026. Patients diagnosed between January 2012 and December 2020 were retrospectively identified, and the study protocol was subsequently reviewed and approved by the Ethics Committee. Because this was a retrospective study involving the evaluation of anonymized data, the Ethics Committee of the Affiliated People's Hospital of Jiangsu University waived the requirement for patient consent.

Results

Patient Characteristics

In the entire cohort of 1231 patients, the median follow-up duration was 120.6 (IQR 95.1-145.5) months, estimated using the reverse Kaplan-Meier method, and the median DFS was 109.6 (IQR 85.0-134.2) months, based on the conventional Kaplan-Meier analysis. In the external TCGA validation cohort, the median follow-up duration was 49 (IQR 24.8-76.1) months, whereas the median DFS was 42.5 (IQR 24.6-73.6) months. At 5 and 8 years, 122 and 70 patients, respectively, remained at risk in the TCGA cohort (Figure 1).

‎

The baseline characteristics of patients in the training, internal validation, and external validation cohorts are summarized in Table 1. Across the 3 cohorts, significant differences were observed in most variables, except laterality, reflecting heterogeneity in patient populations and data sources. Nevertheless, the relatively comparable distributions between the training and internal validation cohorts supported model development and internal validation, while the independent external cohort provided an assessment of model generalizability. Importantly, the model maintained satisfactory predictive performance despite the heterogeneity of the external cohort, which is encouraging for its potential clinical applicability. However, these findings should be interpreted cautiously, and further validation in larger, multicenter cohorts is warranted to confirm the robustness and generalizability of the model.

As shown in Figure S1 in Multimedia Appendix 1, most pairwise associations among variables were weak to moderate (Cramér V<0.30), with stronger correlations limited to clinically related features. Nevertheless, all covariates included in the multivariable analyses had variance inflation factors <3 (Table S2 in Multimedia Appendix 1), indicating no substantial multicollinearity.

Model Development and Performance

Multivariable Cox regression analysis (Table 2) identified several variables significantly associated with DFS. Patients who received endocrine therapy had a lower hazard of DFS events than those who did not (hazard ratio [HR] 0.28, 95% CI 0.18‐0.43; P<.001). Compared with patients with T1 disease, those with T2, T3, and T4 disease had higher hazards of DFS events. Similarly, compared with N0 disease, N1, N2, and N3 disease were associated with higher hazards of DFS events (HRs 1.68, 2.08, and 2.73, respectively). ER status was also significantly associated with DFS in the multivariable model, with the direction of the estimated association differing from that observed in the univariable analysis after adjustment for other covariates, including postoperative endocrine therapy. Kaplan-Meier curves were subsequently plotted for all variables that demonstrated significance in the univariable analysis (Figure S2 in Multimedia Appendix 1).

Table 3 presents the predictive performance of 5 ML models in the development, internal validation, and external validation cohorts, evaluated based on discrimination and overall prediction error. In the development cohort, all models demonstrated good discriminatory performance. The C-indexes of the XGBoost (AFT) and LightGBM models exceeded 0.80. Consistently, their AUC values at 3, 5, and 8 years were also above 0.80 (Figure 2A-C).

‎

In the internal validation cohort, a slight attenuation in predictive performance was observed across all models. Nevertheless, XGBoost (AFT) maintained competitive discrimination (C-index=0.784, 95% CI 0.727‐0.833) and stable time-dependent AUCs at 3, 5, and 8 years (0.828, 0.822, and 0.789, respectively) (Figure 2D-F). In the external validation cohort, the XGBoost (AFT) model consistently demonstrated superior predictive performance, achieving a C-index of 0.750 (95% CI 0.654‐0.843). Furthermore, the model exhibited robust discriminative ability over time, with AUC values of 0.791, 0.762, and 0.768 at 3, 5, and 8 years, respectively (Figure 2G–I). Regarding other performance metrics, the XGBoost (AFT) model demonstrated optimal performance in predictive accuracy (0.740) and clinical utility (sensitivity=0.784, NPV=0.700) (Table S3 in Multimedia Appendix 1).

Calibration curves were generated for the 3 models with the highest C-index values in the external validation cohort (Figure S3 in Multimedia Appendix 1). In the training cohort, the 3-, 5-, and 8-year calibration curves of XGBoost (AFT) closely aligned with the ideal 45 degree line, indicating good agreement between predicted and observed risks. Clinical utility was assessed using decision curve analysis (Figure S4 in Multimedia Appendix 1), which revealed that the predictive curves for all cohorts remained well above the 2 baseline strategies. Although the decision curves intersected, XGBoost (AFT) provided comparable or higher net benefit across a clinically relevant range of threshold probabilities, supporting its potential utility in practice.

To assess the incremental predictive performance of the XGBoost (AFT) ML model relative to the traditional Cox proportional hazards regression model, we calculated the net reclassification improvement (NRI) and integrated discrimination improvement (IDI) at 3, 5, and 8 years of follow-up (Table S4 in Multimedia Appendix 1). For 3-year survival prediction, the XGBoost (AFT) model yielded an NRI of 0.338 (95% CI 0.180‐0.492, P=.01) and an IDI of 0.132 (95% CI 0.074‐0.188, P<.001). At 5 years, the corresponding values were an NRI of 0.311 (95% CI 0.172‐0.458, P=.01) and an IDI of 0.166 (95% CI 0.095‐0.235, P<.001). Similarly, for 8-year survival prediction, the XGBoost (AFT) model demonstrated a significant improvement over the Cox model, with an NRI of 0.265 (95% CI 0.113‐0.452, P=.02) and an IDI of 0.187 (95% CI 0.095‐0.281, P<.001). Across all evaluated time horizons, both NRI and IDI were statistically significant, indicating that the XGBoost (AFT) model achieved superior risk stratification and predictive discrimination compared with the conventional Cox model.

Overall, XGBoost (AFT) achieved the highest C-index among the evaluated models while maintaining competitive performance on other evaluation metrics across the cohorts. Considering its discriminative performance, predictive accuracy, and SHAP-based interpretability, XGBoost (AFT) was selected for subsequent interpretability analysis.

Model Interpretability

To further elucidate the interpretability of the XGBoost (AFT) model, we used SHAP to quantify the contribution of each clinical feature to the model’s predictions of breast cancer DFS.

The SHAP summary plot (Figure 3A) illustrates the global importance of variables ranked according to their mean absolute SHAP values. pT stage emerged as the most influential feature, with a mean absolute SHAP value of 0.460, followed by endocrine therapy, LNR, pN stage, Ki-67, and chemotherapy (Table S5 in Multimedia Appendix 1). These findings indicate that pathological staging, particularly lymph node–related variables, together with treatment-related and molecular features, contributed substantially to the model’s DFS predictions.

‎

As shown in Figure 3B, the SHAP beeswarm plot illustrates the direction and magnitude of feature contributions to the model predictions across different variable levels. Red dots representing N0 status and a low LNR (LNR<0.20) were predominantly distributed on the left, indicating contributions toward lower model-predicted risk, whereas those representing N3 disease and a high LNR (LNR>0.65) were mainly concentrated on the right, indicating contributions toward higher model-predicted risk. Similarly, endocrine therapy was characterized by predominantly leftward SHAP values, indicating a contribution toward lower model-predicted risk within the model. Red dots for T1 tumors were predominantly clustered on the left, whereas those for T3 tumors were mainly distributed on the right, indicating greater contributions toward higher model-predicted risk with increasing tumor stage.

As shown in Figure 3C, the SHAP force plot illustrates individual-level feature contributions to the XGBoost (AFT) model prediction for a randomly selected patient. The baseline model output was 0.078, whereas the final model output for this patient was 0.355, reflecting the combined contributions of the individual features. Thus, the difference between the patient-level prediction and the baseline output represents the aggregate contribution of the included features to the model prediction. Importantly, SHAP values describe feature contributions to model predictions and should not be interpreted as causal effects on DFS.

Overall, SHAP analysis provided global and individual-level explanations of the XGBoost (AFT) model predictions and improved the transparency and interpretability of the model. The observed contribution patterns were also broadly consistent with established clinical knowledge in breast oncology.

Development of DFS Prediction Composite Index

SHAP feature importance analysis of the preliminary XGBoost (AFT) model showed comparable contributions of LNR and pN to model predictions. To further evaluate the robustness of the model to the inclusion of these variables, 2 additional XGBoost (AFT) models were constructed by separately excluding LNR or pN. Neither model showed a significant difference in predictive performance compared with the initial model containing both variables (Table S6 in Multimedia Appendix 1). Given their similar contributions and the comparable predictive performance observed in the sensitivity analyses, pN was retained because of its established role in TNM staging and greater clinical applicability, whereas LNR was excluded to improve model parsimony and facilitate clinical implementation.

The DFS score model demonstrated good discriminative ability across the training, internal validation, and external validation cohorts, with C-indices of 0.735, 0.725, and 0.715, respectively. The time-dependent AUC values at 3, 5, and 8 years consistently exceeded 0.70. In comparison, the traditional TNM staging system had a C-index of 0.618, with time-dependent AUC values of 0.625, 0.593, and 0.614 at 3, 5, and 8 years, respectively (Table 4).

Patients were stratified into low- and high-risk groups using the median risk score as the cutoff. DFS differed significantly between the 2 groups (log-rank test, P<.001), with a higher model-predicted risk corresponding to shorter DFS (HR 5.44, 95% CI 4.04‐7.31; P<.001). A significant trend toward shorter DFS was also observed across increasing risk levels (trend test, P<.001). The 5-year DFS rates were 91.7% and 60.2% in the low- and high-risk groups, respectively. Kaplan-Meier curves showed early separation between the risk groups, and the curves remained separated during follow-up (Figure 4).

‎

Discussion

Principal Findings

Based on a clinical cohort of patients with breast cancer from our institution, this study evaluated 5 ML algorithms, with Cox regression as a conventional statistical reference, to develop a model for postoperative DFS prediction. SHAP analysis was used to characterize the contributions of individual features to model predictions and improve model interpretability. A composite DFS score was further developed for individualized prognostic risk stratification [18-20].

Our model was constructed using a clinical cohort characterized by long-term follow-up, with a median DFS of 109 (IQR 85.0-134.2) months (approximately 9.1 y) among patients. This duration significantly exceeds that of many retrospective studies focusing on 5-year survival rates, highlighting a key advantage of our predictive model. An extensive body of research indicates that 1 to 3 years after surgery represents the peak period for breast cancer recurrence, and the risk of recurrence increases gradually after 5 years, while extending adjuvant therapy beyond this period can improve patient outcomes [21-23]. Moreover, high-risk patients exhibit a higher risk of recurrence across all time intervals, suggesting that the period beyond 5 years after surgery is also a critical phase for disease progression [24,25]. The long-term follow-up of this cohort not only supports the assessment of the ML model’s stability in long-term prediction [26] but also enables its evaluation using data from TCGA. Recent studies have further emphasized the importance of dynamic biomarkers and long-term prognostic assessment in breast cancer, including changes in hormone receptor and HER2 status after neoadjuvant chemotherapy [27]. In the present study, the 5-year outcome was considered the primary time horizon for evaluating model performance because of its clinical relevance and the relatively more mature follow-up available at this time point. The 8-year analysis was additionally performed to explore the potential long-term applicability of the model; however, because of the limited long-term follow-up in the TCGA external validation cohort, the 8-year results should be interpreted as exploratory findings rather than confirmatory evidence.

Our results demonstrated that the XGBoost (AFT) model exhibited the best performance among the 5 ML models in the development, internal validation, and external validation sets. In the development set, the C-index for the XGBoost (AFT) model was 0.808. The AUC values for predicting 3-, 5-, and 8-year DFS were 0.833, 0.840, and 0.824, respectively. Notably, the 5-year AUC was 0.840 in the development cohort and remained 0.762 in the external validation cohort, supporting the model’s ability to provide meaningful medium-term DFS predictions across different populations. Furthermore, this model demonstrated optimal performance in both the external and internal validation sets, indicating its predictive stability. Although the 8-year AUC was also evaluated, this result should be interpreted cautiously because the TCGA external validation cohort had limited long-term follow-up. Therefore, the 8-year findings are considered exploratory and should not be interpreted as definitive evidence of sustained long-term predictive performance. In contrast, the traditional TNM staging model had a C-index of 0.618 in the external validation cohort, with corresponding 3-, 5-, and 8-year AUC values of 0.625, 0.593, and 0.614, respectively. These results indicate that the predictive model established in this study outperformed the traditional TNM staging model for predicting DFS. Compared with the Cox model, the XGBoost (AFT) model demonstrated more pronounced improvements in the NRI and IDI, thereby potentially enhancing its clinical utility. These findings are consistent with recent studies demonstrating the potential of ML and SHAP-based approaches to improve prognostic and treatment-response prediction in breast cancer [28,29].

The traditional TNM staging system is the clinical gold standard for assessing breast cancer prognosis. However, it is limited to anatomical features such as tumor size, lymph node involvement, and distant metastasis. This limitation of anatomical staging overlooks the prognostic impact of clinicopathological factors, including patient age, histology, molecular subtype, and adjuvant therapy, thereby limiting its ability to precisely evaluate the influence of clinical heterogeneity and continuous variables on prognosis [30]. In contrast to the linear and proportional hazards assumptions of traditional TNM staging models, ML algorithms can effectively capture nonlinear relationships between features and prognostic end points, as well as complex interactions among features [31]. This may explain the significantly superior predictive performance of the ML models compared with the traditional model in this study. Previous studies have often focused on constructing and validating single-algorithm models due to limitations such as algorithm selection bias and insufficient generalizability [32]. This study systematically compared the predictive performance of 5 ML algorithms, with Cox regression included as a conventional statistical reference, to identify the model with the best overall performance for postoperative DFS prediction in our cohort. Model development and internal validation were conducted using a real-world clinical cohort from our institution, whereas external validation was performed using the independent TCGA cohort. This validation strategy provided an assessment of the model’s generalizability to an independent patient population and data source. Recent studies have similarly explored ML-based prognostic models using TCGA and Gene Expression Omnibus datasets, highlighting the potential of integrating clinical and molecular information for individualized risk prediction [33-35]. Nevertheless, the external validation in the present study should be considered supportive evidence of model generalizability rather than a substitute for prospective multicenter validation.

Although ML models can achieve strong predictive performance, their complex decision-making processes may be difficult for clinicians to interpret, which can hinder their translation into clinical practice [36]. In this study, SHAP analysis was applied to the optimal XGBoost (AFT) model to provide both global and individual-level interpretations of model predictions. The results identified pT, endocrine therapy, pN, LNR, Ki-67, and chemotherapy as among the most important contributors to the model predictions. Consistent with previous studies, higher pT, LNR, pN, and Ki-67 expression, as well as the absence of endocrine therapy and chemotherapy, contributed to higher model-predicted risk [37,38]. Importantly, LNR and pN showed comparable contributions to model predictions. To further evaluate the robustness of their inclusion, sensitivity analyses were performed by separately excluding LNR or pN from the XGBoost (AFT) model. Neither model showed a significant difference in predictive performance compared with the model containing both variables (Table S6 in Multimedia Appendix 1). Given their comparable predictive performance, pN was retained in the final model because of its established role in TNM staging and broader clinical applicability, whereas LNR was excluded to improve model parsimony.

The SHAP force plot further illustrated how individual features contributed to the prediction for a randomly selected patient relative to the model’s baseline output. This individual-level explanation provides a transparent representation of how the included features collectively shaped the model prediction and may facilitate the interpretation of patient-specific risk estimates. Overall, SHAP analysis provided both global and individual-level explanations of the XGBoost (AFT) model predictions, demonstrating that the model’s decision patterns were broadly consistent with established clinical knowledge in breast oncology while improving the transparency and interpretability of the model. This interpretable framework may support individualized prognostic risk assessment and risk stratification in patients with breast cancer.

Accurate prognostic risk stratification is important for the individualized management and follow-up of patients with breast cancer [39]. However, the TNM staging system has limited ability to capture prognostic heterogeneity among patients with the same stage and does not incorporate multiple clinical, pathological, and treatment-related factors [40]. In this study, a composite DFS risk score was developed by integrating variables identified through SHAP-based feature importance analysis of the optimal ML model with independent prognostic factors identified by multivariable Cox regression. The risk score was calculated as a linear combination of the selected variables, with the corresponding regression coefficients used as variable weights [41]. The resulting DFS score demonstrated better predictive performance than the traditional TNM staging system. Using the median risk score as the cutoff, patients were stratified into low- and high-risk groups, with significantly different DFS between the groups. These findings indicate that the score can effectively distinguish patients according to their model-predicted prognostic risk and may provide a practical approach to postoperative risk stratification. Its potential clinical value lies primarily in identifying patients with different levels of postoperative recurrence risk and supporting individualized surveillance and further clinical assessment. However, because this was a retrospective observational study, the risk score should not be interpreted as a tool for directly recommending treatment intensification or de-escalation. In particular, the treatment variables represented postoperative treatments actually received, and their predictive importance should therefore be interpreted as prognostic information within the model rather than as causal treatment effects.

Recent research has further expanded the application of ML in breast cancer prognosis by incorporating dynamic biomarkers, multiomics data, and multimodal information. Changes in hormone receptor and HER2 status following neoadjuvant chemotherapy have been associated with subsequent oncological outcomes, suggesting that dynamic biomarker assessment may provide prognostic information beyond baseline measurements [27]. Recent studies have similarly suggested that ML-derived prognostic signatures may provide complementary information to conventional clinicopathological staging by integrating multiple dimensions of disease biology. Importantly, emerging studies have begun to extend this approach beyond conventional clinical variables by integrating molecular, multiomics, pathological, and imaging information [42]. For example, multiomics and ML studies have identified potential prognostic biomarkers and biological pathways, including PLA2G4A- and APOO-related molecular characteristics [34,35]. In addition, recent studies have developed prognostic models using multiple ML algorithms and immune- or metabolism-related molecular features [33]. Multimodal approaches integrating clinical information with pathology-derived features, deep learning features, and MRI radiomics have also shown promising performance for survival prediction in breast cancer [43]. These findings suggest that future iterations of our model could potentially incorporate longitudinal blood biomarkers, digital pathology, radiomics, and genomic or transcriptomic information while retaining SHAP-based interpretability. Explainable AI may provide a useful framework for integrating heterogeneous data sources and presenting complex model predictions in a clinically interpretable manner [29]. Other recent real-world studies, including pharmacovigilance analyses of breast cancer therapies, have further highlighted the value of complementary real-world evidence for evaluating treatment-related outcomes and safety [44]. Although these outcomes were beyond the scope of the present DFS model, they underscore the importance of evaluating predictive models within real-world clinical contexts. Methodological studies in other cancers have likewise demonstrated the potential of integrating ML with immune, transcriptomic, and single-cell information [45]; however, these findings should not be directly extrapolated to breast cancer.

Limitations

This study has several limitations. First, the retrospective single-center design and use of the TCGA cohort for external validation may introduce selection and information biases, whereas differences in patient characteristics and treatment patterns may affect model performance. Although external validation provided supportive evidence of generalizability, prospective multicenter validation is warranted. We are also collaborating with other institutions to collect independent clinical data, including blood-based laboratory indicators and additional clinicopathological variables, for future validation. Second, the model mainly incorporated conventional clinicopathological and treatment-related variables and did not include genomic, transcriptomic, multiomics, or radiomic information, which may provide additional prognostic information. Third, the TCGA cohort had limited long-term follow-up. Thus, although 8-year DFS prediction was evaluated, these results should be considered exploratory, and the 5-year results should be regarded as the principal evidence of model performance. Fourth, treatment variables represented postoperative treatments actually received, and their predictive importance should be interpreted as prognostic information within the model rather than as causal treatment effects. Fifth, the institutional cohort spanned the transition from the seventh to the eighth edition of the American Joint Committee on Cancer staging system during 2012 to 2020. Although the anatomic pT and pN categories used in this study did not undergo substantive changes between the 2 editions and were therefore comparable across the study period, minor staging-related heterogeneity cannot be completely excluded. In addition, Ki-67 was assessed by immunohistochemistry in the institutional cohort, and potential interobserver and methodological heterogeneity in Ki-67 assessment may have affected its classification. In the TCGA-BRCA cohort, standardized immunohistochemistry measurements were unavailable; therefore, MKI67 and TP53 mRNA expression were used as transcriptomic surrogate measures for external validation. These transcriptomic measures are not directly equivalent to immunohistochemistry-derived Ki-67 and p53 protein expression, particularly because TP53 transcript abundance does not directly reflect p53 protein accumulation. Although sensitivity analysis excluding TP53 showed essentially unchanged external validation performance and SHAP contribution patterns, this transcriptomic-protein-level discrepancy should be considered when interpreting the external validation results. Finally, although SHAP enhanced model interpretability, prospective clinical evaluation is needed to establish the clinical utility of the model.

Conclusions

In summary, this study developed and externally validated an interpretable ML model for postoperative DFS prediction in patients with breast cancer. The XGBoost (AFT) model demonstrated better predictive performance than conventional TNM staging, whereas SHAP analysis enhanced model interpretability. The derived risk score enabled prognostic risk stratification and provided a quantitative approach to individualized risk assessment and postoperative follow-up management. Further prospective multicenter validation is warranted before broader clinical application.